Chapter 6 Beta diversity

load("data/data.Rdata")

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  dplyr::select_if(~!all(. == 0)) %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  dplyr::select_if(~!all(. == 0)) %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  dplyr::select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, tree = genome_tree)

beta_q1f <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  dplyr::select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, dist = dist)

6.1 Permanova

#Richness
betadisper(beta_q0n$C, sample_metadata$Origin) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00649 0.0064888 0.2231    999  0.637
Residuals 90 2.61727 0.0290808                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
        Feral  Tame
Feral         0.634
Tame  0.63781      

adonis2(beta_q0n$C ~ Origin, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q0n$C))), 
        permutations = 999, 
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q0n$C))) %>% pull(Location)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Origin	1	0.360299	0.01814891	1.663594	0.23
Residual	90	19.492078	0.98185109	NA	NA
Total	91	19.852377	1.00000000	NA	NA

#Neutral diversity
betadisper(beta_q1n$C, sample_metadata$Origin) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01097 0.010969 0.5763    999  0.454
Residuals 90 1.71289 0.019032                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
        Feral  Tame
Feral         0.456
Tame  0.44974      

adonis2(beta_q1n$C ~ Origin, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$C))), 
        permutations = 999, 
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1n$C))) %>% pull(Location)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Origin	1	0.3587603	0.01623493	1.485257	0.287
Residual	90	21.7392877	0.98376507	NA	NA
Total	91	22.0980479	1.00000000	NA	NA

#Phylogenetic diversity
betadisper(beta_q1p$C, sample_metadata$Origin) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01235 0.012355 0.6076    999  0.444
Residuals 90 1.83012 0.020335                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
        Feral  Tame
Feral         0.447
Tame  0.43775      

adonis2(beta_q1p$C ~ Origin, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1p$C))), 
        permutations = 999, 
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1p$C))) %>% pull(Location)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Origin	1	0.1824935	0.02489257	2.297522	0.157
Residual	90	7.1487493	0.97510743	NA	NA
Total	91	7.3312427	1.00000000	NA	NA

#Functional diversity
betadisper(beta_q1f$C, sample_metadata$Origin) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.1982 0.198217 2.7199    999  0.089 .
Residuals 90 6.5589 0.072877                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
        Feral  Tame
Feral         0.099
Tame  0.10259      

adonis2(beta_q1f$C ~ Origin, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1f$C))), 
        permutations = 999, 
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1f$C))) %>% pull(Location)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Origin	1	0.3017094	0.02299117	2.117898	0.244
Residual	90	12.8211288	0.97700883	NA	NA
Total	91	13.1228382	1.00000000	NA	NA

6.2 Richness diversity plot

beta_q0n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(Origin,Location) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = Origin, fill = Origin, shape = as.factor(Location))) +
      scale_color_manual(name="Origin",
          breaks=c("Tame","Feral"),
          values=c("#6A9AC3","#F3B942")) +
      scale_fill_manual(name="Origin",
          breaks=c("Tame","Feral"),
          values=c("#6A9AC350","#F3B94250")) +
    geom_point(size = 4) +
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) +
    labs(shape="Individual")

6.3 Neutral diversity plot

beta_q1n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(Origin,Location) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = Origin, fill = Origin, shape = as.factor(Location))) +
      scale_color_manual(name="Origin",
          breaks=c("Tame","Feral"),
          values=c("#6A9AC3","#F3B942")) +
      scale_fill_manual(name="Origin",
          breaks=c("Tame","Feral"),
          values=c("#6A9AC350","#F3B94250")) +
    geom_point(size = 4) +
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) +
    labs(shape="Individual")